Synchronizability determined by coupling strengths and topology on Complex Networks 
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We investigate in depth the synchronization of coupled oscillators on top of complex networks with different 
degrees of heterogeneity within the context of the Kuramoto model. In a previous paper [Phys. Rev. Lett. 
98, 034101 (2007)], we unveiled how for fixed coupling strengths local patterns of synchronization emerge 
differently in homogeneous and heterogeneous complex networks. Here, we provide more evidence on this 
phenomenon extending the previous work to networks that interpolate between homogeneous and heterogeneous 
topologies. We also present new details on the path towards synchronization for the evolution of clustering in 
the synchronized patterns. Finally, we investigate the synchronization of networks with modular structure and 
conclude that, in these cases, local synchronization is first attained at the most internal level of organization 
of modules, progressively evolving to the outer levels as the coupling constant is increased. The present work 
introduces new parameters that are proved to be useful for the characterization of synchronization phenomena 
in complex networks. 

PACS numbers: 05.45.Xt, 89.75.Fb 



I. INTRODUCTION 

Studies on the emergence of collective and synchronized 
dynamics in large ensembles of coupled units have been car- 
ried out since the beginning of the nineties in different con- 
texts and in a variety of fields, ranging from biology, ecol- 
ogy, and semiconductor lasers, to electronic circuits jl], 0, 0] . 
Collective synchronized dynamics has multiple applications 
in technology, and is a common framework to investigate the 
crucial features in the emergence of critical phenomena in nat- 
ural systems. For instance, it is a relevant issue to fully un- 
derstand some diseases that appear as the result of a sudden 
and undesirable synchronization of a large number of neu- 
ronal units [4]. Recently, synchronization phenomena have 
also been proved to be helpful outside the traditional fields 
where it applies, for instance, in sociology where it can be 
used to study the mechanisms leading to the formation of so- 
cial collective behaviors HHH). 

Among the many models that have been proposed to ad- 
dress synchronization phenomena, one of the most successful 
attempts to understand them is due to Kuramoto [0, Ht], who 
capitalized on previous works by Winfree {§}], and proposed a 
model system of nearly identical weakly coupled limit-cycle 
oscillators. The Kuramoto (KM) mean field case correspond- 
ing to a uniform, all-to-all and sinusoidal coupling is de- 
scribed by the equations of motion, 



e i = w i + -5^sin(fl i -fl i ) (i = l,...,N). (1) 

i=i 

where the factor is incorporated in order to ensure a good 
behavior of the model in the thermodynamic limit, N — > oo, 
u>i stands for the natural frequencies of the oscillators, and 
K is the coupling constant. Moreover, the coherence of the 
population of N oscillators is measured by the complex order 



parameter, 

1 N 

r(t) exp (i0(t)) = - cxp (i0 3 (t)) , (2) 

where the modulus < r(t) < 1 measures the phase co- 
herence of the population and <fi(t) is the average phase. In 
what follows, we will focus on the synchronization of coupled 
oscillators described by the dynamics Eq. (fTJ, because of its 
validity as an approximation for a large number of nonlinear 
equations and its ubiquity in the nonlinear literature IToll . 

The KM approach to synchronization was a great break- 
through for the understanding of the emergence of synchro- 
nization in large populations of oscillators, in particular it 
presents a second-order phase transition from incoherence to 
synchronization, in the order parameter Eq.© for a critical 
value of the coupling constant. However, a large amount of 
real systems do not show a homogeneous pattern of intercon- 
nections among their parts fTTl [Tltl where the original KM 
assumptions apply. 

Many real natural lfl3i [Till , social lfl5ll and technological 
lfl6l [l7l [l8tl systems conform as networks of nodes with con- 
nectivity patterns that diverge considerably from homogene- 
ity, and are usually characterized by a scale-free degree dis- 
tribution, P(k) ~ fc~ 7 (the degree k is the number of con- 
nections of a node). The study of processes taking place on 
top of scale-free networks has led to reconsider classical re- 
sults obtained for regular lattices or random graphs due to 
the radical changes of the system's dynamics when the het- 
erogeneity of the connectivity patterns can not be neglected 
IT^[2^l2ll[2l[23tl24]l. In this case one has to deal with two 
sources of complexity, the nonlinear character of the dynam- 
ics and the complex structures of the substrate, which are usu- 
ally entangled. A contemporary effort to attack this entangled 
problem was due to Watts and Strogatz, that in 1998, trying to 
understand the synchronization of cricket chirps, which show 
a high degree of coordination over long distances as though 
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the insects where "invisibly" connected, end up with a seminal 
paper 12511 about the small-world connectivity property. This 
work was the seed of the modern theory of complex networks 
Ull [T2I1 . Nevertheless, the understanding of the synchroniza- 
tion dynamics in complex networks still remains a challenge. 

In recent years, scientists have addressed the problem of 
synchronization on complex networks capitalizing on the 
Master Stability Function (MSF) formalism rt26ll which al- 
lows to study the stability of the fully synchronized state 
J23, EE |33. EH H2. |33l H3l. The MSF is the result of 
a linear stability analysis for a completely synchronized sys- 
tem. While the MSF approach is useful to get a first insight 
into what is going on in the system as far as the stability of 
the synchronized state is concerned, it tells nothing about how 
synchronization is attained and whether or not the system un- 
der study exhibits a transition similar to the original KM. To 
this end, one must rely on numerical calculations and explore 
the entire phase diagram. Surprisingly, there are only a few 
works that have dealt with the study of the whole synchroniza- 
tion dynamics in specific scenarios i35l [36l [37L l38l l39l l40ll as 
compared with those where the MSF is used, given that the 
onset of synchronization is reacher in its behavioral repertoire 
than the state of complete synchronization. 

In a previous work l4lll . we have shown how, for fixed 
coupling strengths, local patterns of synchronization emerge 
differently in homogeneous and heterogeneous complex net- 
works, driving the process towards a certain degree of global 
synchronization following different paths. In this paper, we 
extend the previous work to different topologies, even those 
with modular structure, and report more results supporting the 
prev ious claim. First, we extend the analysis carried out in 
114111 to networks in which the degree of heterogeneity can be 
tuned between the two limits of random scale-free networks 
and random graphs with a Poisson degree distribution. Sec- 
ond, in order to get further insights about the role of the struc- 
tural properties on the route towards complete synchroniza- 
tion, we study the same dynamics on top of networks with a 
non-random structure at the mesoscopic level, i.e., networks 
with communities. The results support the usefulness of the 
tools developed and highlight the relevance of synchroniza- 
tion phenomena to study in detail the relationship between 
structure and function in complex networks. 



H. KM MODEL ON COMPLEX NETWORKS 

Let us now focus on the paradigmatic Kuramoto model. In 
order to manage with the KM on top of complex topologies 
we reformulate eq. (Q~|) to the form 

dB- 

-±=LO i + J2 A v A ij s H e 3-0i) {i = l,-,N), (3) 
3 

where Ay is the coupling strength between pairs of connected 
oscillators and Aij is the connectivity matrix (Ay = 1 if 
i is linked to j and otherwise). The original Kuramoto 
model introduced above assumed mean-field interactions so 
that Aij = 1, Vi 7^ j (all-to-all) and Ay = K/N, V«, j. 



The first problem when dealing with the KM in complex 
networks is the definition of the dynamics. In the seminal pa- 
per by Kuramoto |7|], eq. (U), the coupling term in the right 
hand side of eq. ® is an intensive magnitude. The depen- 
dence on the number of oscillators N is avoided by choosing 
Aij = K/N. This prescription turns out to be essential for the 
analysis of the system in the thermodynamic limit N — > 00. 
However, choosing Ay = K/N the dynamics of the KM in a 
complex network becomes dependent on N. Therefore, in the 
thermodynamic limit, the coupling term tends to zero except 
for those nodes with a degree that scales with TV l42tl . 

A second prescription consists of taking Ay = K/ki 
(where fcj is the degree of node i) so that Ay is a weighted 
interaction factor that also makes intensive the right hand side 
of Eq. (0). This form has been used to solve the so-called 
paradox of heterogeneity that states that the heterogeneity in 
the degree distribution, which often reduces the average dis- 
tance between nodes, may suppress synchronization in net- 
works of oscillators coupled symmetrically with uniform cou- 
pling strength l34ll . One should consider this result carefully 
because it refers to the stability of the fully synchronized state 
(see below) not to the whole evolution of synchronization in 
the network. More important, the inclusion of weights in the 
interaction strongly affects the original KM dynamics in com- 
plex networks because it imposes a dynamic homogeneity that 
could mask the real topological heter ogeneity of the network. 

Finally, the prescription Ay = K Q9ll43ll44ll . which may 
seem more appropriate, also presents some conceptual prob- 
lems because the sum in the right hand side of eq. (0 could 
eventually diverge in the thermodynamic limit if synchroniza- 
tion is achieved. To our understanding, the most accurate in- 
terpretation of the KM dynamics in complex networks should 
preserve the essential fact of treating the heterogeneity of the 
network independently of the interaction dynamics, and at the 
same time, should remain calculable in the thermodynamic 
limit. Taking into account these factors, the interaction Ay 
in complex networks should be inversely proportional to the 
largest degree of the system Ay = K/k max = A keeping 
in this way the original formulation of the KM valid in the 
thermodynamic limit (in SF networks k max ~ N 1 ^ 1 ^ 1 ^). In 
addition, the same order parameter, eq. (0, can be used to 
describe the coherence of the synchronized state. Since k max 
is constant for a given network, the physical meaning of this 
prescription is a re-scaling of the time units involved in the 
dynamics. Note, however, that for a proper comparison of the 
synchronizability of different complex networks, the global 
and local measures of coherence should be represented ac- 
cording to their respective time scales. Therefore, given two 
complex networks A and B with k max = kA and k max = 
ks respectively, the comparison between observable must be 
done for the same effective coupling Ka/^a = K b /&b = A. 
With this formulation in mind eq. (f3]) reduces to 

^=u i + \^A ij wa.(e i -6 i ) (i = l,...,N), (4) 

3 

independently of the specific topology of the network. This 
allow us to study the dynamics of eq. (0|i over different topolo- 
gies in order to compare the results and properly inspect the 
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interplay between topology and dynamics in what concerns to 
synchronization. 



III. HOMOGENEOUS VERSUS HETEROGENEOUS 
TOPOLOGIES 

Recent results have shed light on the influence of the lo- 
cal interactions' topology on the route to synchronization 
rf3ol E(3]. However, in these studies at least two parameters 
(clustering and average path length) vary along the studied 
family of networks. This paired evolution, although yielding 
an interesting interplay between the two topological param- 
eters, makes it difficult to distinguish what effects were due 
to one or other factors. Here, we would like to address first 
what is the influence of heterogeneity, keeping the number of 
degrees of freedom to a minimum for the comparison to be 
meaningful. The family of networks used in the present sec- 
tion are comparable in their clustering, average distance and 
correlations so that the only difference relies on the degree dis- 
tribution, ranging from a Poissonian type to a scale-free distri- 
bution. Later on in this paper, we will relax these constraints 
and study networks in which the main topological feature is 
given at the mesoscopic scale, i.e., networks with community 
structure. 

Therefore, let us first scrutinize and compare the synchro- 
nization patterns in Erdos-Renyi (ER) and Scale-Free (SF) 
networks. For this purpose we make use of the model pro- 
posed in ll45ll that allows a smooth interpolation between these 
two extremal topologies. Besides, we introduce a new param- 
eter to characterize the synchronization paths to unravel their 
differences. The results reveal that the synchronizability of 
these networks does depend on the coupling between units, 
and hence, that general statements about their synchronizabil- 
ity are eventually misleading. Moreover, we show that even in 
the incoherent solution, r = 0, the system is self-organizing 
towards synchronization. We will analyze in detail how this 
self-organization is attained. 

The first numerical study about the onset of synchroniza- 
tion of Kuramoto oscillators in SF networks l39ll revealed the 
great propense of SF networks to synchronization, which is 
revealed by a non-zero but very small critical value A c 
Besides, it was observed that at the synchronized state, r = 1, 
hubs are extremely robust to perturbations since the recovery 
time of a node as a function of its degree follows a power law 
with exponent —1. However, how do SF networks compare 
with homogeneous networks and what are the roots of the dif- 
ferent behaviors observed? 

We first concentrate on global synchronization for the Ku- 
ramoto model Eq. (0). For this we follow the evolution of 
the order parameter r, Eq. (|2), as A increases, to capture 
the global coherence of the synchronization in networks. We 
will perform this analysis on the family of networks gener- 
ated with the model introduced in ll45ll . This model generates 
a one-parameter family of networks labeled by a G [0, 1]. 
The parameter a measures the degree of heterogeneity of the 
final networks so that a = corresponds to the heterogeneous 
BA network and a = 1 to homogeneous ER graphs. For in- 
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FIG. 1: Global synchronization curves r(A) for different network 
topologies labeled by a (a = corresponds to the BA limit and 
a = 1 to ER graphs). The inset shows the region where the onset 
of synchronization takes place. The network sizes are N = 10 4 



and (k) = 6 (M = 3 
introduced in 14511 



10") and were generated using the model 



termediate values of a one obtains networks that have been 
grown combining both preferential attachment and homoge- 
neous random linking so that each mechanism is chosen with 
probabilities (1 — a) and a, respectively. It is worth stress- 
ing that the growth mechanism preserves the total number of 
links, Ni, and nodes, N, for a proper comparison between dif- 
ferent values of a. Specifically, assuming the final size of the 
network to be N, the network is build up starting from a fully 
connected core of mo nodes and a set 5(0) of N — too un- 
connected nodes. Then, at each time step, a new node (not 
selected before) is chosen from 5(0) and linked to to other 
nodes. Each of the to links is attached with probability a 
to a randomly chosen node (avoiding self- connections) from 
the whole set of N — 1 remaining nodes and with probabil- 
ity (1 — a) following a linear preferential attachment strategy 
114711 . After repeating this process N — too times, networks in- 
terpolating between the limiting cases of ER (a = 1) and SF 
(a = 0) topologies are generated ll45ll . Furthermore, with this 
procedure, the degree of heterogeneity of the grown networks 
varies smoothly between the two limiting cases. 

The curves r(A) for several network topologies ranging 
from ER to SF are shown in FigQ] We have performed exten- 
sive numerical simulations of eq. © for each network sub- 
strate starting from A = and increasing it up to A = 0.4 with 
S\ = 0.02. A large number (at least 500) of different net- 
work realizations and initial conditions were considered for 
every value of A in order to obtain an accurate phase diagram. 
The natural frequencies to*, and the initial values of 0i were 
randomly drawn from a uniform distribution in the interval 
( — 1/2, 1/2) and (— ir, 7r), respectively. 

FigQ] reveals the differences in the critical behavior as a 
function of the substrate heterogeneity. The global coherence 
of the synchronized state, represented by r, shows that the on- 
set of synchronization first occurs for SF networks. As the net- 
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TABLE I: Topological properties of the networks used in this work 
and critical points for the onset of synchronization obtained from a 
FSS analysis (Eq. $5$). The topological quantities reported are the 
result of an average over 1000 network realizations, (k) = 4 and 
TV = 10 4 have been set for all networks. Standard deviation of the 
mean values for A c is ±2 units in the last significant digit. 

{k ) kmax A c 
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work substrate becomes more homogeneous the critical point 
A c shifts to larger values and the system seems to be less syn- 
chronizable. On the other hand, it is also clear that the route to 
complete synchronization, r = 1, is faster for homogeneous 
networks. That is, when A > X c (a) the growth rate of r in- 
creases with a. To inspect in depth the critical parameters of 
the system dynamics we perform a finite size scaling analysis. 
This allows to determine with precision the curve A c (a) and 
study the critical behavior near the synchronization transition. 
We assume a scaling relation of the form 



r = N-"f(NP(\-\ c )), 



(5) 



where f(x) is as usual a universal scaling function bounded as 
x — > ±oo and v and f} are critical exponents to be determined. 
The detailed analysis performed for both SF and ER topolo- 
gies shows that the critical value of the effective coupling, A c , 
corresponds in scale-free networks to Af F = 0.051, and in 
random networks to Xf R = 0.122, accordingly with Fig[TJ In 
both cases, the transition strongly recalls the classical transi- 
tion of the original KM 10] with a critical exponent near 1/2 
for the SF network ll39tl . For intermediate values of a, the re- 
sults show that the critical point shifts to larger values as the 
degree of heterogeneity increases. They are shown in Table U 
together with some topological properties of the networks. 

The differences between ER and SF topologies observed 
when looking at global patterns of synchronization motivate 
a more detailed study of the synchronization onset for both 
topologies. The original work by Kuramoto pointed out that 
at the onset of synchronization small clusters of locked oscil- 
lators emerge and that the recruitment of more oscillators into 
these clusters as the coupling is increased makes it larger the 
global coherence r of the system. Obviously the emergence 
of these clusters would depend on the underlying topology 
which drives the possible configurations that locked oscilla- 
tors would eventually form. To see how this initial coherence 
is achieved we propose a new order parameter, r^nfe- This 
parameter measures the local construction of the synchroniza- 
tion patterns ll48ll and allows for the exploration of how global 
synchronization is attained. We define 
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FIG. 2: Evolution of the control parameters r and runk as a func- 
tion of the coupling strength for networks generated with the model 
introduced in 145I1 . corresponding to a = 0.0, 0.25, 0.5, 0.75 and 
1.0. The size of the networks is TV = 10 3 and their average degree 
is (k) = 6. The exponent of the SF networks increases from 7 = 3 
(a = 0). 



being T; the set of neighbors of node i. The parameter ru n k 
measures the fraction of all possible links that are synchro- 
nized in the network. The averaging time At should be taken 
large enough in order to obtain good measures of the degree of 
coherence between each pair of physically connected nodes. 
Besides, r^nfe is computed after the system relaxes at some 
large time t r . Note that in the limit of all-to-all coupling the 
information provided by ru n k is exactly the same that the one 
provided by r because in this case ru n k oc r 2 . Therefore, no 
additional information would be provided by this new param- 
eter in the all-to-all case. Here, however, it turns out to be the 
key parameter to characterize how synchronization emerges at 
a local scale. 

In Fig[2]we represent the evolution of both order parame- 
ters, r and Tymk, as a function of the coupling strength A for 
several values of a. The behavior of Hmfc shows a change in 
synchronizability between ER and SF and provides additional 
information to that reported by r. Interestingly, the nonzero 
values of ninfc f° r A < A c indicate the existence of some local 
synchronization patterns even in the regime of global inco- 
herence (r w 0). Right at the onset of synchronization for 
the SF network limit, its runk value deviates from that of the 
ER recovering the known result about the synchronization of 
SF networks for lower values of the coupling. In this region, 
while the synchronization patterns continue to grow for the 
ER network at the same rate, the formation of locally syn- 
chronized structures occurs at a faster rate in the SF network. 
Finally, when the incoherent solution in the ER network desta- 
bilizes, the growing in its synchronization pattern increases 
drastically up to values of ru n k comparable to those obtained 
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in SF networks and even higher. For intermediate values of 
a, the results show that the effect of varying the heterogene- 
ity of the underlying network is twofold. On one hand, the 
more heterogeneous the network is, the smaller the values of 
A needed for the onset of synchronization. Conversely, the in- 
crease in the degree of heterogeneity results in larger values of 
A in order to achieve complete synchronization. In short, as 
the heterogeneity is increased, the onset of synchronization is 
anticipated, but at the same time, the appearance of the fully 
synchronized state is delayed. 

These results undoubtedly point out that statements about 
synchronizability are dependent on the coupling strength 
value. To shed new light on this phenomenon, we have stud- 
ied the characteristics of the synchronization patterns along 
the evolution of runk- Following the usual picture, synchro- 
nization patterns are formed by pairs of oscillators, physically 
connected, whose phase difference in the stationary state tends 
to zero. In order to determine which pairs of nodes are truly 
synchronized we should determine the coherence of their dy- 
namics. Note that eq.© is the average dynamical coherence 
between every pair of linked nodes and then the synchroniza- 
tion degree of every pair of connected oscillators can be writ- 
ten in terms of a symmetric matrix 
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(7) 



Then one has to analyze each matrix term in order to la- 
bel a link as synchronized or not. As introduced above, 
from the computation of m n k one determines the fraction of 
physical links that are synchronized so that one would expect 
that 2ru n k ■ Ni elements of the matrix T> are D^ = 1, while 
the remaining elements are D^ = 0. However, this is not the 
real situation since the network dynamics is not well defined 
in terms of a fully synchronized cluster and a set of completely 
incoherent oscillators. On the other hand, the worst scenario 
would be found if there were 27Y; elements of matrix V so 
that Dij = rimk, implying that all the physically connected 
pairs are equally synchronized and hence the parameter ru n k 
could not be interpreted as the fraction of links that are dy- 
namically coherent and no information about the topological 
patterns of synchronization could be extracted from matrix T>. 
The situation found is not as simple as the former possibility 
and not so dramatic as the latter. The contributions of 
the Ni elements of matrix T> that correspond to physical links 
can be ordered from the highest to the lowest one. We have 
checked that for two situations, corresponding to the onset of 
synchronization (A = 0.05) and when high global coherence 
(A = 0.13) is observed for a SF network, synchronized links 
can be clearly identified. For the onset of synchronization, a 
subset of nearly 20% of links displaying coherent dynamics 
with high degree of synchronization, XV, > 0.8, is well sep- 
arated from the behavior of the remaining links as a dramatic 
decrease of Dij takes place. In this sense, it is clear that the 
dynamics of a 20% of the possible pairs can be regarded as 
synchronized which is in agreement with the obtained value 
TUnk = 0.25 for A = 0.05 and support that although macro- 
scopic coherence is not observed (r ~ at this point) the sys- 
tem is seen to walk towards it. For A = 0.13 (ru n k — 0.82) 
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FIG. 3: Evolution of the number of synchronized clusters N c and the 
synchronized giant component size GC as a function of ru n k f° r the 
the different topologies considered. Small values of ru n k correspond 
to values of A for which r ~ 0. Despite r being vanishing and 
hence no global synchronization is yet attained, a significant number 
of clusters show up. This indicates that for any A > the system 
self-organizes towards macroscopic synchronization. The network 
parameters are as in Fig. [2] 



a plateau of nearly 75% of links is observed, thus revealing 
the high degree of global coherence, r ~ 0.7, at this point. 
Therefore, the shape of the ranked curves confirm that 
Tank gives the fraction of synchronized links and thus the lat- 
ter allows to obtain information about synchronized patterns 
from T>. 

To determine exactly which pairs of nodes are regarded as 
synchronized, the matrix T> is filtered using a threshold T such 
that the fraction of synchronized pairs equals ru n k- In this 
way T is a moving threshold so that if T>ij > T oscillators i 
and j are considered synchronized. The value of T depends 
on the particular realization and is determined by means of 
an iterative scheme starting from T = 1, Decreasing it with 
ST = 0.01 one computes the amount of links that fulfills the 
condition. In this way, the value of T progressively decreases 
and more pairs of oscillators are chosen. The process lasts 
until T is such that the fraction of chosen links is equal to 
the desired value ru n k previously computed from V. Finally, 
when the synchronized links are identified the clusters of syn- 
chronized nodes are reconstructed. 

Figure [3] represents the number of synchronized clusters 
and the size of the giant component (GC) as a function of 
funk for the same values of a used in Fig. [2] The local in- 
formation extracted from it points to a novel feature of the 
synchronization process that is not possible to derive from 
Figs[T]and|2j and that is unexpected. The emergence of clus- 
ters of synchronized pairs of oscillators (links) in the networks 
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FIG. 4: Giant synchronized components for several values of A in the two limiting cases of the different topologies studied (ER and SF). The 
size of the underlying networks is small (N — 100 nodes), in order to have a sizeable picture of the system. Note that for the SF case links and 
nodes are incorporated together to the GC, while for the ER network, what is added are links between nodes already belonging to the GC. 



shows that for values of A for which the incoherent solution 
r = is stable, the networks have developed a largest clus- 
ter of synchronized pairs of oscillators involving 50% of the 
nodes, and an equal number of smaller synchronization clus- 
ters. From this point on, the behavior of both GC and N c 
depends on the specific value of a. When heterogeneity dom- 
inates, the GC grows and the number of smaller clusters goes 
down, whereas for less heterogeneous networks the growth 
of GC is more abrupt and nodes are incorporated to it more 
faster. Moreover, the results highlight the fact that although 
heterogeneous networks exhibit more coherence in terms of 
r and ru n k, the microscopic evolution of the synchronization 
patterns is faster in homogeneous networks, being these net- 
works far more locally synchronizable than the heterogeneous 
ones once A > A c . 

The observed differences in the behavior at a local scale are 
rooted in the growth of the GC. For homogeneous topologies, 
many small clusters of synchronized pairs of oscillators (note 
in Fig|3]the large number of clusters formed when a 15% of 
the links are synchronized) merge together to form a GC when 
the effective coupling is increased. This coalescence of many 
small clusters results in a giant component made up of almost 
the size of the system once the incoherent state destabilizes. 
On the other hand, for heterogeneous graphs, the growth of the 
giant component is more smooth and the oscillators form new 
pairs starting from a core made up of half the nodes of the net- 
work. That is, in one case (ER-like networks), almost all the 
nodes of the network takes part of the giant component from 
the beginning and latter on, when A is increased, what is added 
to the GC are the links among these nodes that were missing 
in the original cluster of synchronized nodes. For SF-like net- 
works, the mechanism is the opposite. Nodes are added to the 
GC together with most of their links, resulting in a growth of 
funk much slower than for the homogeneous topologies. 

The above picture is confirmed in Figj4] where we have rep- 
resented the evolution of the local synchronization patterns of 
the giant components in ER and SF networks for several val- 
ues of A l49ll . It is clear that when r ~ the two networks fol- 
low different paths toward synchronization. In particular, the 
giant component for the SF network seems to retain the topo- 
logical features of the substrate network, while this is not the 
case for the ER network (for instance, the small-world prop- 



erty is clearly lacking). 

This study about the patterns of self-organization towards 
synchronization reveals that the quantitative difference about 
the macroscopic behavior, shown by the computation of the 
evolution of the global coherence r for ER and SF networks, 
has its roots on a qualitatively different route at the micro- 
scopic level of description. The use of the new parameter r^nfe 
which involves the computation of the degree of coherence be- 
tween each pair of linked nodes is a useful tool for describing 
such differences. Moreover, the results suggest that the degree 
of heterogeneity of the network is the key ingredient to explain 
the two different routes observed. The technique developed to 
extract the synchronization patterns allows the analysis of the 
topological features of such clusters of nodes. We can com- 
pute the average measures of relevant quantities such as the 
clustering coefficient or the degree distribution, and see how 
these magnitudes evolve from the uncoupled limit, where no 
synchronization occurs, to the coherent regime where the syn- 
chronized network coincides with the underlying substrate. It 
is then relevant to explore the regions where the onset of syn- 
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FIG. 5: Evolution of the ratio between the clustering coefficient of 
the giant synchronized cluster, (c sync }, and that of the substrate net- 
work (cnetwork), as a function of A for the two limiting cases of BA 
and ER networks. Network parameters are those used in figure|2] 
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chronization takes place and characterize topologically these 
emergent synchronized clusters. 

In Fig|5] the evolution of the average clustering coeffi- 
cient (csync) of the giant synchronized cluster referred to 
(cnetwork) m the underlying network, is plotted as a function 
of A for both the BA and ER networks. It is worth mentioning 
that the results depicted in the figure have been computed tak- 
ing into account that nodes with degree 1 does not contribute 
to the clustering coefficient of the GC, as c is not properly 
defined for these nodes. The results are illustrative of the lo- 
cal organization of synchronized nodes. The figure shows that 
for both topologies the clustering decreases as the coupling 
is increased beyond their respective A c or, in other words, as 
the giant component grows by the addition of new synchro- 
nized pairs of nodes. However, the effects of the two different 
routes to complete synchronization observed for ER and SF 
networks are well appreciated from the results. For the hetero- 
geneous network there is a smooth decrease of the clustering 
coefficient for A > Af F and the effects of the emergence of 
global coherence are not dramatic in what refers to the behav- 
ior of (c sync }. This is because in this case the giant component 
mainly grows by recruiting new synchronized nodes and their 
links. On the other hand, for the ER graph the behavior ob- 
served for A < Xf R , i.e. when no macroscopic coherence 
is observed, is interrupted by a sudden jump near its critical 
value. In fact, for A > Xf R the clustering of the synchro- 
nized cluster quickly approaches the value of (c) of the sub- 
strate network. This effect becomes clear if one has in mind 
the coalescence of small clusters, which happens around the 
critical point for ER graphs. In fact, taking into account the 
giant synchronized component on ER for A < Xf R , implies 
to consider one of the several disjoint synchronized clusters of 
similar sizes that are in this region. Moreover, the coalescence 
process leads to the formation of a giant cluster that contains 
almost all the nodes of the network (see Fig.O, but a number 
of links significantly smaller. Hence, when the clusters col- 
lapse into a much larger one, the topological features change 
dramatically as observed from the evolution of the clustering 
coefficient. 

All the results reported above point out that the ultimate 
reason behind the two different routes to complete synchro- 
nization is the heterogeneous character of the SF network and 
the role played by the hubs. The natural cohesion that hubs 
provide to SF network prevents the existence of independent 
macroscopic clusters of synchrony as occurs for ER networks. 
It is then interesting to study how these hubs participate in the 
formation of the final synchronized state. For this, we first 
study the evolution with A of the composition of the synchro- 
nized cluster in terms of the degree of its components. In l4lll . 
we reported the probability that a node with degree k belongs 
to the giant synchronized cluster as a function of its degree 
k and the coupling A for the SF network. This probability 
turns out to be an increasing function of k for every value 
of A and therefore the more connected a node is, the more 
likely it takes part in the cluster of synchronized links. In 
particular, the results confirm the hypothesis made above that 
the hubs participate from the very beginning on the formation 
of the synchronized cluster. A similar result was obtained in 




FIG. 6: (color online) The plot shows the fraction of links that a node 
with degree k belonging to the synchronized cluster shares with other 
nodes of the same synchronized cluster. This fraction ki nt /k is plot- 
ted as a function of k and A. The figure shows how the hubs progres- 
sively incorporate their neighbors to the synchronized component as 
A grows. The network is SF with parameters as those used in Fig[2] 
and a = 0. 

l5lll . where Zhou and Kurths studied the hierarchical organi- 
zation in complex networks, using the MSF and a mean-field 
approach in the weak coupling limit. 

The above characterization of the synchronized cluster in 
terms of the degree of its component can be completed study- 
ing their effective degree, ki n t . The effective degree of a syn- 
chronized node is the number of links it shares with other 
nodes belonging to the same synchronized cluster. Obviously, 
at the complete synchronized regime a node with degree k will 
have ki„t = k. We have plotted in Fig|6]the quantity kint/k 
(the fraction of links that a node has with synchronized neigh- 
bors) as a function of A and the degree k of the nodes (a = 0). 
The results reveal that although hubs are the first to take part 
of the synchronized cluster, their neighbors are progressively 
incorporated to the cluster as A grows. Besides, if a node with 
small k is synchronized the probability that its neighbors are 
also synchronized grows very fast with A which is an effect 
of the network topology. These results further support the 
statement about the essential role played by the hubs in the 
recruitment of oscillators into the synchronized group and in 
the emergence of complete synchronization in SF networks. 



IV. SYNCHRONIZATION IN STRUCTURED NETWORKS 

In light of the results of the above section we have extended 
the study beyond unstructured networks to structured or mod- 
ular networks. This is a limiting situation in which the lo- 
cal structure may greatly affect the dynamics, irrespective of 
whether or not we deal with homogeneous or heterogeneous 
networks and then they constitute a perfect framework for test- 
ing the new order parameter runk introduced in the last sec- 
tion. 

Many complex networks in nature are modular, i.e. com- 
posed of certain subgraphs with differentiated internal and ex- 
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ternal connectivity that form communities lfl2i l50ll . The use 
of modular networks where a proper comparison in synchro- 
nizability can be performed (same number of nodes and links) 
restricts us to the consideration of synthetic structured net- 
works. To this end, we make use of a common benchmark 
of random network with community structure, first proposed 
by Newman ll52tl considering one hierarchical level and later 
extended to several hierarchical levels 1371 l38ll . 

The modular network structure we build is as follows: in a 
set of N nodes, we prescribe n compartments that will rep- 
resent our first community organizational level, and m com- 
partments, each one embedding four different compartments 
of the first level, that define the second organizational level of 
the network. The internal degree of nodes at the first level Zj ni 
and the internal degree of nodes at the second level z i7l2 keep 
an average degree zi ni + Zi, l2 + z out = (k) so that these net- 
works are strictly homogeneous in the sense of the degree dis- 
tribution, P(fc) = 6(k— (k)). Networks with two hierarchical 
levels are labeled as - Zm 2 , e.g. a network with i-j means 
i links with the nodes of its first hierarchical community level 
(more internal), j links with the rest of communities that form 
the second hierarchical level (more external) and ((h) —i — j) 
links with any community of the rest of the network. 

Synchronization processes on top of modular networks of 
this type have been recently studied as a mechanism for com- 
munity detection ll37il53ll . In ll37ll . the authors studied the sit- 
uation in which starting from a set of homogeneous (in terms 
of the natural frequencies) Kuramoto oscillators with differ- 
ent initial conditions the system evolves after a transient of 
time to the synchronized state. It was shown that the commu- 
nity structure is progressively unveiled at the same time the 
system's dynamics evolves toward the coherent state and the 
synchronization is attained. In particular, the nodes belonging 
to the first community level are the first to get synchronized, 
subsequently the second level nodes achieve the frequency en- 
trainment and finally the whole system shows global synchro- 
nization. 

Here we adopt a different perspective since we will con- 
sider as previously a set of non-identical Kuramoto oscilla- 
tors with random assignment of natural frequencies and hence 
the final degree of system's synchronization will depend on 
the strength of the coupling. It is then interesting to study 
how the degree of synchronization evolves as a function of 
A and whether the coherence between nodes is progressively 
distributed following the hierarchy imposed by the underly- 
ing topology. For this, we make use of the order parameters 
r, eq. (fJJ, and ru n k, eq ©, to characterize the synchroniza- 
tion transition on two slightly different modular networks with 
two well defined hierarchical levels, 13 — 4 and 15 — 2, being 
this difference the cohesion of the internal community core, 
13 links out of 15 possible neighbors or 15 links (i.e., all-to- 
all) at the most internal level. Both networks have N = 256 
and (k) = 18. FigJT] shows the results for both kinds of net- 
works revealing that the path towards synchronization as a 
function of the interaction is again affected by the structure. 
They also show that the information provided by ru n k is es- 
sential to unveil the synchronization process. While the global 
synchronization parameter r is reflecting that the 13 — 4 struc- 
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FIG. 7: (color online) Evolution of r (top) and ru n u (bottom) as a 
function of A for structured modular networks. The networks are syn- 
thetically built with an a priori community structure. The network 
size is 256 nodes and the number of links is 4608. We prescribe 
16 compartments that will represent our first community organiza- 
tional level, and four compartments each one embedding four differ- 
ent compartments of the above first level, that define the second or- 
ganizational level of the network. Each node has 18 links distributed 
between its first community level, the second, and the whole network 
at random. The network 13-4 has 13 internal connections in its first 
hierarchical level, 4 external connections in its second hierarchical 
level, and 1 connection with any other community in the network. 
The generation of the 15-2 structure is equivalent. The curves show 
that although 13-4 has always a better global synchronization, 15-2 
has better local synchronization as shown by r^nfc- 



ture globally synchronizes always better, ru„k tells us again 
about the local synchronization. It shows that local synchro- 
nization is indeed favored in the 15 — 2 structure since ru n k 
is larger for this topology for small values of A where the sys- 
tem is locally forming synchronized clusters. This result, not 
captured by the macroscopic indicator r, is expected since the 
internal cohesion of communities at the first hierarchical level 
is larger for the 15-2 than for the 13-4. The evolution of rn n k 
shows that when the coupling A is increased the number of 
links synchronized in the 13 — 4 network becomes larger than 
in the 15 — 2 structure revealing that complete synchroniza- 
tion is then favored by the presence of more external links 
connecting the first level communities. 

Fig[8] shows the size of the giant component of synchro- 
nized clusters and the number of them as a function of A. An 
interesting effect of the community structure of the networks 
and of the dynamics of the synchronization process is revealed 
in the figure. Right at the value of A where the onset of global 
coherence takes place, the size of GC suddenly falls to in- 
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FIG. 8: (color online) Size of the largest synchronized cluster GC 
(top) and number of clusters N c (bottom) for the same networks of 
Fig-El S ee the text for details. 
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FIG. 9: Number of links (Nu n k s ) that connect synchronized nodes 
in each of the three levels in the community hierarchy of the net- 
works (1 means inner layer). The numbers are normalized by the 
total number of links at each level in each network. 



crease again at larger values of the coupling strength. Addi- 
tionally, note that this point coincides with that correspond- 
ing to a change in the concavity of the ru n k{\) curves. This 
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FIG. 10: (color online) We represent the degree of synchronization 
between pairs of connected nodes for several values of the coupling 
A in a 13 — 4 modular network (with two organizational levels) of 
N = 256 nodes. The color code denotes the value of the aver- 
aged (over different initial conditions) filtered matrix (23y) 6 [0, 1]. 
The values of the coupling are (from left to right and top to bottom) 
A = 0.011, 0.026, 0.032, 0.035, 0.038, 0.046, 0.210 (corresponding 
to full synchronization). The pictures show that the order of synchro- 
nization is given by the organizational levels. The first community 
level is the first one to get synchronized, subsequently, second level 
nodes attain synchronization for a larger value of A and finally the 
full synchronized state is reached when outer links have (Dij) — 1, 



change at the microscopic level is due to the readjustment of 
links that connect synchronized nodes. In fact, as Fig.|9]illus- 
trates for both networks, in this region of A values, the number 
of links connecting synchronized nodes of the third level de- 
creases while the number of those ascribed to the second level 
raises. That is, the synchronization process takes place in such 
a way that the first to synchronize are the nodes of the inner 
community level, then the second and so on until the whole 
network gets synchronized. The relevant fact is that in order 
for ru n k and r to grow, the nodes and links of the second level 
adjust their phases at the expense of those of the outer layer, 
the third level. This is also reflected in the number of clusters 
of synchronized links (N c ), i.e., the network appears like if 
the nodes of the third level were "temporarily" disconnected. 
Moreover, as the 13 — 4 network has more links connecting the 
first and the second hierarchical levels, Nn n k S 2 raises faster in 
this network than in the 15 — 2. 

We have further inspected the synchronization path in mod- 
ular networks. This can be easily done and visualized by the 
representation of the filtered matrix T>. It implies to reassign 
the values of matrix T> so that T>ij = 1 if 2?y > T , and 
T>ij = otherwise. Plotting this filtered matrix for different 
values of the coupling A one can easily determine which links 
are the first to synchronize since the form of the adjacency ma- 
trix (that includes all the physical links between nodes) is also 
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easy to interpret because of its nested structure. FigfTUlshows 
how the community structure determines the internal organi- 
zation of the system in the route towards full synchronization 
for the 13 — 4 network. For this study we have computed the 
value of the filtered matrix T> for a number of initial condi- 
tions and then took its average value so that (T>ij) e [0, 1] 
accounts for the synchronization strength of the network link 
The results point out that link synchronization depends 
on the organizational level they belong to. Those connect- 
ing nodes belonging to the same first level community are the 
fastest (in terms of the coupling strength A) to reach full syn- 
chronization. For larger values of A full synchronization is at- 
tained progressively for the subsequent organizational levels. 
Then, one can conclude that the inner the link is the faster it 
gets synchronized in agreement with previous studies reported 
above JH. 

V. CONCLUSIONS 

In this paper we have explored several issues about syn- 
chronization in complex networks of Kuramoto phase oscil- 
lators. Our main concern has been the study of the synchro- 
nization patterns that emerge as the coupling between non- 
identical oscillators increases. We have described the degree 
of synchronization between each pair of connected oscillators. 
The use of a new parameter, ru n k, allows to reconstruct the 
synchronization clusters from the dynamical data. We have 
studied how the underlying topology (ranging from homo- 
geneous to heterogeneous structures) affects the evolution of 
synchronization patterns. The results reveal that the route to- 
wards full synchronization depends strongly on whether one 
deals with homogeneous or heterogenous topologies. In par- 
ticular, it has been shown that a giant cluster of synchroniza- 
tion in heterogeneous networks comes out from a unique core 
formed by highly connected nodes (hubs) whereas for homo- 
geneous networks several synchronization clusters of similar 
size can coexist. In the latter case, a coalescence of these clus- 
ters is observed in the synchronization path which is macro- 
scopically manifested by the sudden growth of global coher- 
ence. Another important effect of the underlying topology 
is manifested in an anticipated onset of global coherence for 
heterogeneous networks with respect to more homogeneous 
topologies. However, the latter reaches the state of full syn- 
chronization at lower values of the coupling strength, there- 
fore showing that statements about synchronizability of com- 
plex networks are relative to the region of the phase diagram 
where they operate. Additionally, we have shown that these 
systems are seen to organize towards synchronization even 
when no macroscopic signs of global coherence is observed. 



Finally, the framework of structured networks has provided 
a useful benchmark for testing the validity of the new param- 
eter runk and the information obtained from the computation 
of matrix T>. The results obtained by means of these quan- 
tities allow to conclude that for modular networks synchro- 
nization is first locally attained at the most internal level of 
organization and, as the coupling is increased, it progressively 
evolves toward outer shells of the network. The latter pro- 
cess is however achieved at the expense of partially readjust- 
ing some pairs of synchronized nodes between the inner and 
outer community levels. Besides, we have obtained evidences 
that a high cohesion at the first level communities produce a 
high degree of local synchronization although it delays the ap- 
pearance of the global coherent state. 

This study has extended the previous findings about the 
paths towards synchronization in complex networks fiUl . and 
provides a deeper understanding of phase synchronization 
phenomena on top of complex topologies. In general, the 
work supports the idea that in the absence of analytical tools 
to confront the resolution of non-linear dynamical models in 
complex networks, the introduction of new parameters to de- 
scribe the statistical properties of the emergence of local pat- 
terns is needed as they give novel and useful information that 
might guide our comprehension of these phenomena. On 
more general grounds, this work adds to other recent find- 
ings |p4l l55ll about the topology emerging from dynamical 
processes. The evidences that are being accumulated point to 
a dynamical organization, both at the local and global scales, 
that is driven by the underlying topology. Whether or not this 
intriguing regularity has something to do with the ubiquity 
of complex heterogeneous networks in Nature is not clear yet. 
More works in this direction are needed, but we think that they 
may ultimately lead to uncover important universal relations 
between the structure and function of complex natural systems 
that form networks. Another issue to explore in future works 
concerns the behavior of non-linear dynamical systems on top 
of directed networks ll56ll . which will allow deeper insights 
into the behavior of natural systems. 
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